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ABSTRACT 

We present an update to the multiphase SPH galaxy formation code by Scannapieco et 
al. We include a more elaborate treatment of the production of metals, cooling rates 
based on individual element abundances, and a scheme for the turbulent diffusion 
of metals. Our SN feedback model now transfers energy to the ISM in kinetic and 
thermal form, and we include a prescription for the effects of radiation pressure from 
massive young stars on the ISM. We calibrate our new code on the well studied 
Aquarius haloes and then use it to simulate a sample of 16 galaxies with halo masses 
between 1 x lO'^^ and 3 x 10^"^ Mq. In general, the stellar masses of the sample agree 
well with the stellar mass to halo mass relation inferred from abundance matching 
techniques for redshifts z — — A. There is however a tendency to overproduce stars 
at z > 4 and to underproduce them at z < 0.5 in the least massive haloes. Overly 
high SFRs at z < 1 for the most massive haloes are likely connected to the lack of 
AGN feedback in our model. The simulated sample also shows reasonable agreement 
with observed star formation rates, sizes, gas fractions and gas-phase metallicities at 
z = — 3. Remaining discrepancies can be connected to deviations from predictions 
for star formation histories from abundance matching. At z = 0, the model galaxies 
show realistic morphologies, stellar surface density profiles, circular velocity curves and 
stellar metallicities, but overly flat metallicity gradients. 15 out of 16 of our galaxies 
contain disk components with kinematic disk fraction ranging between 15 and 65 %. 
The disk fraction depends on the time of the last destructive merger or misaligned 
infall event. Considering the remaining shortcomings of our simulations we conclude 
that even higher kinematic disk fractions may be possible for ACDM haloes with quiet 
merger histories, such as the Aquarius haloes. 
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1 INTRODUCTION 

In the standard paradigm of cosmic structure forma- 
tion, galaxies form through coolin g and condensation 
of gas within dark m atter haloes (|White fc ReesI 119781 : 
iFall fc EfstathiouH igsdl. CoUisionless A''-body simulations of 
the dark matter component have been able to reproduce the 
observed large-scale structure of the universe with high ac- 
curacy in the A Cold Dark Ma tter (ACDM) version of this 
paradigm (jSpringel et al.l 120061 ). Semi-analytic models rely- 
ing on these simulations and simple analytical prescriptions 
for the baryonic component are capable of repro ducing the 
detailed systematics of the galaxy population (jGuo et al.l 
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I2OIH ). To properly understand the complex dynamical in- 
teractions between gas, stars and dark matter in galaxies 
requires however cosmological hydrodynamical simulations. 

The complexity of properly modeling the many 
baryonic astrophysical processes which play a role in 
the formation of galaxies has led to a long list of mod- 
els over the last two decades (e.g. iNavarro fc White! 



I199J; _lNavarro & Steinmotz) 
Gov ernato et al. 2007; Scann apieco et a l 



es ( e.g. 
19971: ' 



2003: 



Abadi et al 
2009; Agertz et al 



2011). Despite significant recent progress (e. g. I Sales et al 
,201a ICuedes et all I2OIII : I Brook et al] |2012| ). these sim- 
ulations continue to be plagued by a range of problems, 
with dilTerent codes often producing very different galaxies 
for the same initial conditions (Ok a moto et al, 2003 : 
iKeres et al.ll201ll : IScannapieco et aLll2012l ). 



M. Aumer et al. 



The formation of realistic disc galaxies has been shown 
to be especially problematic. The cooling and condensation 
of too much low-angular momentum gas (overcooling) and 
the loss of angular momentum from baryons to the dark 
matter have been a problem since t he first efforts to simu - 
late cosmological galaxy formation (|Navarro &i Benall99lh . 
Once disks have formed they have been shown to be suscep- 
tible to destruction by major mergers (IToomre fc Toomre 
1972 ). by the infall of satellite galaxies IIToth fc Ostriker 



1992 ) and by accretion of misaligned gas (|Scannapieco et al 



200E )(CS09 hereafter). Moreover, the r eorientation of disks 
I Aumer fc White! l2013l : lOkamotol 120131) a nd disk instabili- 



ties (e.g. lNoguchilll999l : lAumer et al.ll2010l ) can enhance the 
bulge-to-disk ratio. 

In addition to these problems, stellar mass to halo 

mass relations from abundance matching techniques have 

shown that the majority of simulatio ns over-produce stars 

Guo et al.ll2O10l: ISawala et al.ll201ll ). especially at high z 



(CS09, iMoster et al.l l2013) 



Modeling the injection of energy fr om supernova ex- 
plosions into the surroun ding gas (e.g. IScannapieco et al.l 
120061 : IStinson et al.l 120061 ) has been shown to be a possi- 
ble mechanism to prevent overly efficient gas cooling, for 
driving large-scale outflows, for removing low angular mo- 
mentum mat erial and thus for producing more realistic disk 
galaxies fe.g. lScannapieco et al.ll2008l : lBrook et al.ll201ll ). In- 
terestingly, simulations applying rather weak feedback have 
been sho wn to be more successful in reproducing early-type 
galax ies IjNaab et al.ll2007l : IOser et al.ll2O10l : I Johansson et al 



2012). Apart from sup ernovae, AGN (e. g. ISpringel et al 



2005 ) and cosmic rays IjUhlig et al.l |2012| ) have been con- 



sidered as possible sources of feedback. 

Simulations including empirical models of momentum- 
driven winds have been shown to improve the agree- 
ment with observed propertie s of the galaxy population 
and the intergalactic med ium (JQppenheimer fc Davg 120061 : 
lOppenheimer et al.l I2OI0I ) . Recently, the input of momen- 
tum and energy from massive young stars in form of stel- 
lar winds and radiation pressure prio r to their explo s ion as 
SNe has been studied in more detail. iHopkins et al.l (|201ll ) 
showed that momentum input into the ISM from radia- 
tion pressure may play a key role in regulating star forma- 
tion. IStinson et al.l (120131 ) included thermal feedback from 
young stars and thus were able to achieve significantly bet- 
ter agreemen t of st ar formation histories with observations. 
lAgertz et al.l (|2012l ) presented the most complete model so 
far of mass, momentum and energy feedback from stars dur- 
ing all stages of their evolution, and concluded that radia- 
tion pressure from young stars has the strongest effect on 
the ISM. 

Apart from feedback, a more realistic modeling of 
star formation and the ISM has been shown to help in 
making simulated gala xies with more realistic properties. 
iGovernato et al.l (|201Cll ) found that increasing the thresh- 
old density used in the modeling of star formation to more 
realistic values leads to more concentrated star formation, 
more efficient feedback energy input into the ISM, and the 
formation of galax ie s with higher dis k fract ion (see also 
iGuedes et all I2OIII ). IChristensen et all l|2012l ) have shown 
that the implementatio n of a m odel for the formation of 
molecular hydrogen (see lGnedin et al. 2009 ) can amplify the 
effect of higher threshold densities. 



In lAumer fc Whitd l|2013l) (AW13 hereafter), we have 
shown that it is possible to form a realistic disk galaxy within 
a tria:xial, substructured and growing CDM halo, as long as 
the disk forms predominantly at late times {z < 1), the an- 
gular momentum vector of the infalling gas has a roughly 
constant and appropriate orientation and there are no ma- 
jor mergers or other significant changes in the configura- 
tion of the halo. In the Aquila Code Comparison Project 
(|Scannapieco et al.ll2012l . CS12 hereafter) one of the haloes 
considered in AW13 was simulated with 15 different galaxy 
formation codes yielding 15 model galaxies with widely vary- 
ing properties. Although several models compared well to a 
subset of the observations considered, none of the models 
yielded a fully realistic disk galaxy. In this paper, we intend 
to show how modifications to the models of CS09, motivated 
by the findings of CS12 and some of the recently discussed 
solutions to problems discussed above, can lead to the for- 
mation of significantly more realistic disk galaxies. We also 
use our model to study a sam ple of haloes, which has been 
previously studied by lOser et al . (2010) in the context of 
the formation of massive, early-type galaxies. The combined 
sample comprises haloes with very quiet merger histories as 
well as low-2 mergers, and is thus well suited to be compared 
to a range of objects in the observed galaxy population. 

Our paper is organized as follows: In Section 2 we de- 
scribe the updates we applied to our code. In Section 3 we 
introduce the sample of initial conditions we use. In Sec- 
tion 4 we discuss the star formation histories of our galax- 
ies. In Section 5 we analyze their kinematical and structural 
properties. In Section 6 we compare our sample to observed 
scaling relations. Finally, we conclude in Section 7. 



2 THE CODE 

For our simulations, w e use the T REES P H code GADGET- 
3, last des cribed in ISpringell (|2005l ). IScannapieco et al.l 
(|2005l . 12009 ) introduced models for stellar metal production, 

metal line cooling, star formation, SN feedback and a mul- 
tiphase gas treatment to be used with GADGET-3, which 
we use as the basis for our code update. 



2.1 Multiphase model and star formation 

We begin with a description of the parts of the model, which 
have remained unchanged. The code is unique in decoupling 
SPH particles with very different thermodynamic properties 
by preventing them from being SPH neighbours. Particle i 
decouples from particle j if 



Ai > 50. X Ai 



and 



/iij < Cij. (1) 

Here A is the entropic function (jSpringel fc Hernauistl2002l ). 
Cij is the pair-averaged sound-speed and ^lij is the relative 
velo city of the partic l es pro jected onto their separation vec- 
tor. iMarri fc White! (|2003! ) showed that the second condi- 
tion is needed to avoid decoupling in shock waves, which 
can lead to unphysical effects. This multiphase gas treat- 
ment has been implemented to make a realistic co-existence 
of hot and cold phase gas (as observed in the ISM) possible 
in SPH. It has also been shown to allow a more realistic 
modeling of energy and metal injection from stars to the 
distinct components of the ISM. 
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To model the formation of stars, we assume that gas 
particles are eligible for this process, if their density is 
above a th reshold density nth. W hereas CS09 used nth ~ 
0.03 cm~^. lGovernato et al.l ([20101) have argued that signif- 
icantly higher threshold densities are needed to form real- 
istic galaxies. While we can confirm their conclusions, we 
caution that the effect of varying this parameter depends 
significantly on the details of the applied feedback and ISM 
model. For the model and the resolution applied in this 
work, we use a value of nth ^ 3 cm~^, fo r whic h we find 
the best results. As argued bv Guedes et al.l (|201ll ). a signif- 
icantly lower value leads to less efficient feedback and higher 
bulge fraction. However, a significantly higher threshold in 
our model leads to the formation of bound stellar clumps, 
which can sink to the centre because of dynamical friction 
and thus also enhances the bulge fraction. Our value lies be- 
tween the value of 1 cm~^ applied for the Gaso line model for 
halo A q-C-5 in CS12 and 5 cm~^ applied bv lGuedes et al.l 
l|201ll) in a higher resolution simulation. We note that these 
values are still two orders of magnitude below the density 
of molecular clouds, however our resolution is too coarse to 
model these. 

For particles with n > nth and an overdensity p/p > 
2000 (where p is the cosmic mean density) which lie in a 
convergent flow, a star formation rate of 



dp* 

dt 



n 



Pgaa 

L-dyn 



(2) 



is assumed. Here stellar and gas densities are represented by 
p* and Pgas and tdyn = 1/ ^jAirGpg^s is the local dynamical 
time for the gas particle. We choose a star formation effi- 
ciency r] = 0.04, in the range of values typically used (see 
e.g. various models in CS12). The typical timescale for star 
formation is thus t^t = l/rj tdyn- Star particles are created 
stochastically with one gas particle being turned into one 
star particle of the same mass. 



2.2 Metal Production and Cooling 

To account for metals, we explicitly trace the mass in the 
elements H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe for all gas 
and star particles. Our model includes chemical enrichment 
from SNH, SNIa and AGB stars. Each star particle repre - 
sents a stellar population characterized by a lKroupal (120011 ) 
IMF with lower and upper mass limits of 0.1 and 100 Mq. 
We assume that stars more massive than 8 Mq ex- 
plode as SNH. All SNH are modeled in one event at an 
age r = 3 Myr as the energy release from SNH peaks at 
this time and the most massive stars are supposed to have 
the strongest effect on the ISM. For the calculation of the 
mass returned to the ISM in the variou s elements we use 
the m etal-dependent yields provided by IChieffi fc Limongil 
(| 2004 1. The uncertainties in yields and thus in the pre- 
dictions of s imulations are s i gnific ant, as was for example 



dictions ot s imulations are s i gmnc ant, as was tor example 
discussed by Wiersma et al.l ((20091). It has been argued by 
IPortinari et al . (1998) that adjustments by factors of 2 for 
certain elements can help improve the agreement with ob- 
servations. Indeed, we find that halving the iron yield follow- 
ing their findings leads to a qualitatively better agreement 
of metallicities as discussed in the following sections, which 
is why we apply their suggested corrections. We note that 
apart from this, we have not studied the variation of our 



results with different yield sets, IMFs etc. to optimize our 
results. 

For the element pr oduction by SNIa, w e assume the 
model W7 presented bv llwamoto et al.l (|l999r ). We apply a 
delay time distribution, which declines with age r of a stel- 
lar p opulation as r~^, as proposed bv iMaoz fc Mannuccil 
( |2012i ). We also adopt their suggested normalization of 2 
SNIa per 1000 Mq of stars formed and that the first SNIa 
explode at r = 50 Myr. The corresponding masses of ele- 
ments are returned to the ISM in 50 Myr time-steps in our 
model. 

To account for the mass recycling in the winds of asymp- 
to tic giant branch stars, we use the metal-dependent yields 
of'Karaka^ (|2010h . Together with the assumed IMF and us- 
ing lifetimes dependent on stellar mass and metallicity, we 
can thus calculate the mass of the considered elements re- 
leased during the time intervals considered for the SNIa en- 
richment. 

Chemical elements are distributed to the gaseous neigh- 
bours of the star particles, where neighbours are weighted 
according to their distance from the star particle using an 
SPH kernel. To account for our multiphase treatment of the 
gas component, the returned (metal) mass is split between 
10 hot and 10 cold neighbours. We give 50 per cent to the 
hot and cold phase each for all three different types of metal 
production sites. We have tested making this fraction de- 
pendent on the age of the stellar particle, but found that 
observed abundance ratios are best reproduced by our sim- 
ple choice. For this purpose the cold phase gas is defined by 
r < 8 X 10** K and n > 4. x 10"^ cm'^. Note, that we have 
reduced this density limi t by a factor of ~ 100 compared to 
IScannapieco et al.l (|2006r ) as we found that a higher value 
can have a destructive effect on extended, low-density gas 
disks due to energetic feedback (see below). 

The metallicities of the SPH particles are used to cal- 
culate the c ooling rates of t he ga s. We apply the rates 
presented by IWiersma et al.l (|2009l ) for optically thin gas 
in ionization equilibrium. These rates are calculated on an 
element-by-element basis and take into account the effects 
of photo-ionizatio n from a uniform redshi ft-dependent ion- 
izing background (jHaardt fc Madaul 120011 ). The rates thus 
depend on redshift, gas density, temperature and chemical 
composition. 



2.3 Metal Diffusion 

In many standard im plementations of chemica l enrichment 
in SPH (also true for IScannapieco et al.ll2005l ). the metal- 
licity of a particle can only change by enrichment. This can 
lead to situations, where gas particles with similar thermo- 
dynamic properties, but very different metallicities live next 
to each other. This occurs e.g. when a galactic wind parti- 
cle travels through unenriched IGM. IWiersma et al.l (|2009r ) 
suggested that smoothed metallicities should be used to get 
rid of these situations. Although this improves the model- 
ing by smoothing out differences, it also leads to ffuctuating 
metallicities for individual particles for the IGM situation 
discussed above. 

In the ISM, once metals have been released from stars, 
turbulent motions of gas are responsible for their spread- 
ing. Including a corresponding model for th e turbu l ent di f- 
fusion of metals was already suggested bv iGroomI (|l997h . 
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but until recently, most galaxy formation SPH codes did 
not consi der this process. Recent i mplementations were pre - 
sent ed inlMartmez-Serrano et all (120081 ) . iGreif et all (|2009l ') 
and lShenet all (|2010l ). 

The diffusion equation for a metal concentration c 
(metal mass per total mass) of a fluid element with den- 
sity p and a diffusion coefficient D is 



dc 
di 



-V ■ (DVc) 



(3) 



where d/dt is th e Lagrangian derivative. 
IClearv fc MonaghanI (| 19981 ') gave the SPH formulation 
of the diffusion equation as 



(4) 



At 


= "^Kij [a ~Cj), 

3 


whei 


e 

TUj 4DiDj Vij ■ ViWij 


iiij 


PiPj{Di + Dj) r?. 



(5) 



Here quantities with subscripts i and j correspond to neigh- 
bour particles, m is the particle mass, Wtj is the SPH kernel 
and Vi^ is the separati on vector with absolute rij. 

Grcif ct alj (|2009|) argued for the use of an integrated 
equation assuming that the change in c is small over Ai: 



(6) 



(7) 



As we want to conserve the total metal mass, we modify the 
equation for a pairwise exchange of metals. For the metal 
mass p,i — CiTUi we get 



Ci{to + At) = 


-- c.ito)e^^' + ^il-c 


With 


and B = N KijCj 



A/i, = J2 ^''^ = H [2 



1 



-e'^^') -TK^oicJ-c)) 



A 



,(8) 



where the factor 1/2 was included to account for the fact 
that most pairs of neighbours are considered twice and fiij 
is correspondingly subtracted from particle j. To avoid de- 
pendence on the ordering of particles all changes Apij are 
calculated for all pairs of neighbours before the metal masses 
of all particles are updated. This procedure is applied at ev- 
ery time-step for all active particles using the standard SPH 
neighbour searches. A neighbour particle j can be inactive, 
so that the corresponding pair of particles only appears once. 
Should the particles still be neighbours at the next active 
time-step of j, the larger Ai compensates for that. Clearly 
this formalism includes a num ber of appro x imati ons, but ap- 
plied to tests as discussed in iGreif et al.l (|2009r ) , we find a 
similar accuracy. 

This leave s the d etermination of the diffusion coefficient 
D. iGreif et al.l (|2009|) argued for Di = 2piaihi, where a is 
the velocity dispersion of gas particles within its smoothing 
kernel characterized by the smoothing length /t,. IShen et al.l 
l|2010l ') argued that Di = O.OSpi \Ski\ hf based on the trace- 
free tensor Ski (for details see their paper) is a better choice 
as it yields no diffusion for purely rotating or compressive 
flows. We have tested both ideas and found that for a fixed 
test setup the main differen c e is th e strength of the diffusion 
coefficient with iGreif et al.l (|2009l ) predicting values higher 



by a factor ~ 20. For cosmological simulations, we find 
that the Shen et al. configuration yields better r esults when 
compa ring to observations. As was noticed by IShen et al.l 
(|201CI ), diffusion leads to outfiowing particles losing met- 
als to the circumgalactic medium and subsequently also to 
higher gaseous and stellar metallicities in the galaxy. For the 
Greif et al. value for D this effect is much stronger than for 
the Shen et al. value and makes galaxies lie above the mass- 
metallicity relation. However when we use Di = Q.lpiOihi 
this criterion loses significance. For consistency, we use D as 
suggested by Shen et al. in the simulations presented in this 
paper. 



2.4 Thermal and Kinetic feedback 



As in IScannapieco et al.l (|200a) we assume that each SN 
ejects an energy of lO^^erg into the surrounding ISM. As 
for the metals we split this energy in halves and give those 
to the 10 nearest hot and cold gas neighbour particles as 
defined above. However, we now split the energy between a 
kinetic and a thermal part (see also Agortz et al. 2012). 

As discussed above we know how much mass Am is 
transferred from a star particle to a neighbour gas parti- 
cle at a given feedback time-step. We also know how many 
SN events are represented and thus the total energy that 
is available to be transferred A_Etot. We assume that the 
momentum of the outfiowing material Ap = AmVout , where 
Uout is the velocity of the outflowing material, which we set 
to be 3000 kms"'^ and pointing radially away from the star 
towards the gas particle, is conserved and transferred to the 
gas particle. This leads to a change in kinetic energy of the 
gas particle A_Ekin. 

The remaining energy is considered to be thermalized, 
A_Ethcr m = Ai?tot — A_Ekin. For thls energy we follow the 
ideas of IScannapieco et al.l (|200fJ). The fraction transferred 
to a hot particle is instantly added to its thermal energy. 
Instead, for cold particle we accumulate the energy from SNe 
events in a reservoir. Only when the accumulated energy is 
high enough, so that it can become a hot particle, the energy 
is released ('promotion'). To define 'hot' in this context, we 
search for neighbour particles with entropies high enough 
to be decoupled from the cold particle in our multiphase 
scheme and calculate their mean entropy Ahot . If a particle 
has less than 5 such neighbours within 10 smoothing lengths, 
we set ^hot ~ j4th, where j4th correspon ds to T = 1.6 x lO^K 
and n = 2 X 10"^cm"^ ('seeding', see ISawala et al.ll201ll '). 
Moreover, ^hot = ^th is also assumed if ^hot < ^th, so that 
Ath acts as a minimum promotion entropy. 

As A_Ethcrm > AiJkin for the Considered situation, the 
thermal feedback is only mildly weakened. The kinetic feed- 
back however helps breaking up dense clumps of gas and 
thus lowers the cooling rates and the star formation effi- 
ciency in a star-forming disk. Gas fractions increase and 
so does the efficiency of the thermal feedback. Momentum 
is also transferred, when mass is released from AGB stars, 
however masses and velocities are much smaller than in SN 
explosions, which is why this effect is negligible. 



2.5 Radiation Pressure 

iHopkins et al.l l|201ll ) used high-resolution simulations to 
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test the idea that the inclusion of feedback from radiation 
pressure of young massive stars has a comparable or possi- 
bly even stronger effect on t he ISM than SN fee dback (for 
the underlying ideas see e.g. iMurrav et al.lboosl 'l. They de- 
veloped a model for high resolution simulations, in which 
they identify star forming regions and can thus use their 
properties to calculate the effect of radiation pressur e as a 
function of these local properties. IStinson et al.l l|2013l ) have 
introduced a model, where instead of a momentum trans- 
fer, a thermal energy transfer is modeled, assuming that 10 
per cent of the radiated energy from massive young stars 
is thermalized in the surrounding ISM. lAgertz et alj l|2012l ) 
presented a model for radiation pressure, which assumes that 
each star particle represents a sample of star forming clus- 
ters. The effect is however then dependent on their choice 
of the maximum cluster mass. 

We parametrize the rate of momentum deposition to 
the gas as 



c 



(9) 



where tir is the infrared optical depth and L{t) is the UV- 
luminosity of the stellar population a nd we have set add i- 
tional efficiency parameters to 1 (see lAgertz et al.l 120121 ). 
This equation states that all UV photons are scattered or 
absorbed by dust, which subsequently re-radiates the en- 
ergy in the infra-red. The IR photons are then scattered 
multiple times before leaving the star forming cloud. We 
construct L{t) by usin g the stellar evolution models for mas- 
sive st ars bylEkstrom et al.l (|2012l ) and assuming a Kroupa 
(|200ll ) IMF. 

More problem a tic is the estimation of tir. If we fol- 
low iHopkins et al.l l|201ll ). tir = Edump^iR with kir ~ 
5 cm^g~^, we effectively have to know the surface density of 
the star forming gas clump. Due to our coarse resolution, we 
do not resolve such clumps. We therefore use the following 
model: 



-TIR : 



Tofip,Z)g{a). 



(10) 



Each gas particle is characterized by its density and metal- 
licity through f{p,Z), whereas the environment of the star 
particle is characterized by g{(j). f{p,Z) basically models 
the dependence on the dust surface density represented by 
a particle: 



fiP,Z) = 



ph 



pthht 



Z 



(11) 



Here we use ph with the smoothing length h as a, measure for 
the surface density represented by a particle and limit the 
effect of this factor at the star formation density threshold. 
We also assume a linear dependence of dust on metallicity 
Z. If we set g = 1, we gain insight into how our model works. 
If we choose tq = 25, which is a value in the range of values 
10-100 found bv lHopkins et al.l l|201ll ) in their models of star 
forming high-z disks, we find that indeed radiation pressure 
significantly reduces high-z star formation rates and brings 
them into better agreement wi th abundance match ing re- 
sults, similar to the findings of IStinson et al.l (|2013l ). How- 
ever, the low-2 star formation efficiencies become too low 
and the final galaxy masses are too small indicating that 
radiation pressure is too strong at later times (see Figure [l] 
and discussion in Section [4. II). 
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0664-4X 


213.74 


3.62 


7.37 


LO 


0616-4X 


235.77 


3.62 


7.37 


LO 



Table 1. Overview over the haloes studied in this paper, includ- 
ing the name of the halo, its virial mass Mvir, the dark matter 
particle mass m^^ and the initial gas particle mass rrigas. Under 
' Origin' we distin guish between haloes from 'CS' (CS09) and 'LO' 
llOser et al.ll2010h . 



However low-z disks show significantly lower gas ve- 
locity dispersions than gas-rich, turbulent high-z disks 
(JGenzel et al.ll2008l l. For a typical Milky- Way type halo, we 
find the average gas velocity dispersion cr is ^ 40 kms^ 
at 2 ~ 2, which is in agreement with observed galaxies at 
that time (IForster S chrcibcr ct al. 2009), and at low redshift 
(J ~ 15 — 20 kms~^, which is too high, but explained by our 
coarse resolution. Star forming regions, and thus the values 
of riR, are thus much larger at high z. If we assume that 
the Jeans mass is a valid estimator for clump mass, then 
the mass scales as Md oc a"^ (we ignore the dependence on 
density as we are not modeling densities above the star for- 
mation threshold properly) . As is pointed out in the re view 
of observed star forming clouds in lAgertz et al.l l|2012r ) , the 
final radius of clusters depends weakly on mass, which leads 
to our crude estimation 



pH - (- 



(12) 



To estimate a in the environment of the star particle, we de- 
termine the velocity dispersion crkcrnoi of gas particles within 
the smoothing kernel for each gas particle. Then we average 
over the values akcmoi of all the neighbour gas particles to a 
young star particle to avoid strong fluctuations. 

We acknowledge that this line of estimation is not par- 
ticularly stringent, but it qualitatively assures that the effect 
of radiation pressure is stronger in galaxies with higher gas 
velocity dispersions and thus larger star forming regions. 

In our simulations, we model the momentum change 
Ap — pipAt where At is the time-step of the star particle, 
as a continuous force acting on the 10 nearest cold neighbour 
particles during the first 30 Myr in the life of a star particle. 
We use To = 25 for ao — 40 kms~^, so that the effect at 
high-2 reduces star formation rates and is weak at low-z. 
We also limit (/((t) at a value of 4 to avoid overly strong 
forces. By construction, we thus find values of tir ~ 20 at 
~ 2 with extreme values going up to 100, whereas we find 
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significantly lower values of the order 1 — 5 at low redshift. 
We discuss the choice of parameters in Section 14.11 



3 THE SAMPLE 

As initial conditions for our cosmological zoom-in sim- 
ulations we use haloes from the Aquarius Project 
ijSpringel et al.ll2008l . CS09), one of which was also used for 
the Aquila Project (CS 1 2). In addtion, we use a selection of 
haloes from lOser et al.1 (|201C| ). 

The Aquarius haloes are a suite of high resolution zoom- 
in resim ulations of regions chosen fro m the Millenium II sim- 
ulation (iBovlan-Kolchin et al.ll2009l ) which follows a cosmo- 
logical box of a side-length 137 Mpc. They were simulated 
from z = 127 assuming a ACDM universe with the following 
parameters: flA = 0.75, Qm = 0.25,r2i, = 0.040, ag = 0.9 and 
Ho = 73 kms~^Mpc~^. The haloes were selected to have a 
similar mass to that inferred for the Milky- Way dark halo 
and to have no neighbour exceedi ng half of their mass w ithin 
1.4 Mpc. For details we refer to ISpringel et al.l l|2008l ') and 
CS09. Because of the second criterion, the haloes have a rel- 
atively quiet low-z merger history. They are thus expected 
to host galaxies with high disk fractions at z — 0. Disks were 
indeed found to form in these haloes by CS09 and CS12, but 
with low disk-to-bulge ratios and several other properties 
that do not compare well to observatio ns. 

The haloes from lOser et al.l ((20101) were selected from 
a simulation of a cosmological box with a side-length of 
100 Mpc, and were simulated from z = 43 assuming 
a ACDM universe with the following parameters: Qa = 
0.74, Qm = 0.26, ^6 = 0.044, as = 0.77 and Ho = 
72 kms~^Mpc~^. Unlike the Aquarius haloes, these haloes 
were not chosen specifically as candidates for disk galaxies 
and the sample thus also contains 2^0 mergers. For de- 
tails we refer to lOser et al.l (|2010h . who resimulated lower 
resolution versions of these haloes. 

The haloes we selected from these two samples are 
listed in Table [T] Our combined sample comprises 16 haloes 
with z — re-simulated virial masses Mvir ranging from 
1.7 X 1O"M0 to 2.4 X 10^^ Mq. The resolutions of the sim- 
ulations are characterized by the initial gas particle masses, 
which lie between 2.87 and 7.37 x 10^ Mq and dark mat- 
ter particle masses between 1.50 and 7.37 x 10® Mq. For 
our re-simulations we applied comoving gravitational soft- 
ening lengths of /ibar ~ 300 pc for gas and stars, and of 
hdm ~ 650 pc for dark matter. We thus use shorter soften- 
ing lengths than in CS09. 



4 STAR FORMATION HISTORIES 

We begin to analyze the outcome of our simulations 
by considering the star formation hi story (SFH) of the 
model galaxies. As has been shown bv lMoster et al.l l|2013l l 
(MNW13 hereafter), the majority of simulations predicts 
stellar galaxy masses that are significantly higher than abun- 
dance matching results suggest at all redshifts 2 = — 4 with 
overproductio n being stronger at high redshift ( for ex cep- 
tions see e.g. IStinson et al.l 120131 : iKannan et al.ll2013l V As 
real disk galaxies form most of their stars at low-z, when 
destructive mergers are rare, it is crucial for simulations to 
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1 \ I 

CS09, 700 pc 
,, , no RP, no kin FB, 700 pc 




RP final model, 300 pc 
RP no Sigma, 300 pc 
RP on SF gas, 700 pc 
RP on SF gas, 300 pc _ 
RP simple, 300 pc 
Mosteretal. (2013) 



4.0 3.0 



2.0 



1.0 



0.5 



0.0 



Figure 1. The evolution of tfie stellar galaxy mass M^^^n^^ vs. 
redshift z in halo AqC for various code configurations (coloured 
lines). Overplotted are the predictions of MNW13 in two ways: 
a) The predicted mass evolution of a typical halo with the z = 
mass of AqC (black dashed line) and b) the predicted evolution 
taking into account the actual halo mass of AqC at six redshifts 
(points with 1 a error-bars). 



produce reasonable SFHs in order to get galaxy structural 
properties right. 

4.1 The effect of changing feedback models on the 
SFH of AqC 

Halo AqC has been studied in CS12 with a variety of dif- 
ferent simulation codes. Although some of these codes pro- 
duce realistic z = masses, none produced a good match to 
abundance matching results for the SFH at all z = — 4. 
Compared to haloes of similar z = mass, AqC is among the 
on es which assemble their dark mass earliest (see Figure Al 
in IScannapieco et al.ll2011h . When developing our feedback 
models, we found that AqC is particularly sensitive to the 
choice of model details, which we demonstrate in this sec- 
tion. We caution, that the effects we discuss here can vary 
from halo to halo and not all conclusions drawn from AqC 
are true for all haloes. Moreover, abundance matching gives 
statistical properties for a galaxy sample and AqC could well 
be an outlier. As we will show below, our calibration method 
is justified by the fact that the feedback model that works 
best for AqC produces a sample of galaxies with reasonable 
properties. 

In Figure [l] we plot the evolution of the stellar mass 
Afstoiiar from 2 = 4 — for various runs of the same initial 
conditions. We compare it to the evolution predicted by the 
abundance matching results of MNW13, who provide two 
different ways of comparison. On the one hand, fitting for- 
mulae for typical SF and accretion histories are presented 
as a function of 2 = halo mass M2oo- We represent these 
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Figure 2. Stellar galactic mass Mgtciiar plotted against halo virial 
mass M200 at redshifts 2 = 0, 1,2 and 3.5. The coloured points 
are the models discussed in this paper. Each colour represents a 
simulation. A colour can appear more than once at high-2 repre- 
senting progenitor haloes. This colour coding is used for figures 
throughout the paper. The gray points are models of the same set 
of haloes run with the code used in CS09. The black lines are the 
abundance matching results of MNW13 with \a regions indicated 
by the dashed lines. 



predictions by the dashed line in Figure [T] On the other 
hand, we know the halo mass at each redshift 2 and can 
thus make use of the fitting formulae for Matciiar/MgooC^)- 
We represent these by six data-points with corresponding 
error-bars from z = 4 to 0. This method yields significantly 
higher predictions for Mstdiar at high z, which refiects the 
early assembly of halo AqC. Clearly the latter predictions 
should be used as a guideline for model calibration. 

The original model by CS09 produces a stellar mass at 
z = 4 that is about an order of magnitude too high com- 
pared to the MNW13 value. It remains more than \a above 
predictions at all redshifts. If we add our updated metal 
production and metal cooling rates (no RP, no kin FB) the 
changes are negligible. Changing to our new thermal and 
kinetic feedback [no RP) also does not improve the SFH at 
high redshifts (the run was stopped at z = 2). We note that 
the change from the CS09 thermal feedback to our new ther- 
mal and kinetic feedback can produce significant changes in 
star formation his tories in other haloes. 

IStinson et al.l ( |2013: ) argued that adding feedback from 
young stars before thei r explosion as SNe wo uld significantly 
shift SF to later times. iKannan et al.l (|2013l ) tested this idea 
on a simulation of a cosmological volume and found good 
agreement with the results of MNW13 at high z. We tested 
various models of feedback from radiation pressure (RP) 
First we consider a simple model for RP {simple RP) where 
all affected gas particles are treated equally independent of 
density, metallicity and gas velocity dispersion (assuming 
an effective optical depth tq = 25, see Equation 10). For 
this model we achieve a reduction of the early SFR so that 
the z = A stellar mass matches predictions. However, the 
feedback is too strong later on and SFRs are drastically re- 
duced so that the z = Q stellar mass is more than 2a low. 
When modifying our model so that the force is only acting 
on star-forming particles with p > pth {RP on SF gas), the 
strength of feedback is reduced especially at high 2, so that 
the SFR is high at 2 > 4, but the agreement at 2 = is 
significantly better, however with low-2 SFRs being still too 
low. A model which includes a metallicity and density de- 
pendence to account for the dust surface density of each gas 
particle as described in Section[23](-RP no sigma) lowers the 
effect of RP at low z when metallicities are higher but the 
densities of the affected gas particles are on average lower, 
as the systems have lower gas fractions. Our full model in- 
cluding the dependence on the gas velocity dispersion a {RP 
final model) models the dependence of RP on the size of star 
forming regions. As a is higher in high-2 galaxies, the effect 
of RP is strengthened at high and weakened at low 2 leading 
to a significant change in the shape of the M* — 2 curve and 
bringing it within la of the predictions of MNW13 at all 
2 < 4; 

iKrumholz fc ThompsonI l|2013l ) have argued that the ef- 
ficiency of the coupling between the radiation field and the 
du sty gas is signifi c antly lower than assumed in the approach 
of iHopkins et al.l (|201 j ). which is the basis of our model 
here. This would imply that radiation pressure is less effi- 
cient in reducing the problem of overly high SFRs in sim- 
ulated galaxies at high 2. We thus caution, that the uncer- 
tainties in the modeled processes remain high and we take 
the agreement of simulated SFH with the results of MNW13 
as a justification for the model we apply. 

We also use Figure [T] to show a dependence of our feed- 
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Figure 3. Star formation rates SFR as a function of redshift z for all our models. The haloes are split into M200 < 5 X U)^^ Mq (left), 
5 X IO^^Mq < M200 < 1-55 X IO^^Mq (middle) and M200 > 1-55 x 10^'^Mq (right). Overplotted are predicted SFRs for typical haloes 
given by MNW13 for the mean, lowest and highest halo mass per panel. 



back on the gravitational softening length e. Going from 700 
to 300 pc baryonic comoving softening length significantly 
reduces the high- 2: SFRs, as shown by two versions of our 
RP on SF gas model. This seems to originate in a higher 
feedback efficiency due to denser structures in the very early 
stages of galaxy formation. Because of this effect, we apply 
Ebaryons = 300 pc for our siiuulations, a significantly lower 
value than in CS09. 

We conclude that the inclusion of a parametrized model 
for feedback from young stars in the from of radiation pres- 
sure is indeed an efficient way to bring star formation his- 
tories into better agreement with results from abundance 
matching techniques. We caution that details of the model- 
ing have a significant effect on the outcome and that these 
effects vary from halo to halo. It is crucial to test the model 
on a set of haloes, as we describe in the next subsection. 



4.2 Applying the model to all haloes 

We now apply our updated galaxy formation code to the full 
sample of haloes described in Section |3] and examine if the 
resulting SFHs agree with observations of the real galaxy 
population. 

We start by discussing stellar masses Matdiar as a func- 
tion of their halo masses M200 at redshifts z = 3.5, 2, 1 and 
in Figure [2] Mstdiar is defined as the mass of stars within a 
radius rgai defined by visually analyzing the spherical stellar 
mass profile Mateiiar(< r) and determining where its radial 
growth saturates. The coloured points in this Figure are for 
the models with our updated code, whereas the gray points 
are for a sample of simulations run on the same set of ini- 
tial conditions with the code version and model parameters 
applied in CS09. Note that the colour coding according to 
halo mass that we apply here is used for all appropriate fig- 
ures throughout the paper. For the haloes 2283 and 0977, 
which are undergoing major mergers at 2 = 0, we use the 
combined mass of the galaxies about to merge, as their dark 
haloes have already merged. At z = 0, the population as a 



whole agrees nicely with the halo occupation modeling by 
MNW13. However, there are trends in the sample that dis- 
agree. The four most massive haloes are 1 — 2a high, whereas 
the lower mass galaxies show a tendency towards overly low 
mass. 

The sample simulated with the CS09 code shows z = 
stellar masses that are generally too high, but in reasonable 
agreement with MNW13 at lower halo mass. The strong 
discrepancy between the old code version and abundance 
matching results becomes better visible at higher redshifts, 
as in Figure [l] The old code version yields stellar masses 
that for all haloes are more than la high at z ^ 1, as was 
also shown for a large number of other codes in MNW13. 

For the models simulated with our new code version the 
situation at 2 ^ 1 is significantly different. Although there 
is still a tendency towards too high mass at 2 ^ 2, all val- 
ues are within la uncertainties of the abundance matching 
results. At 2 = 1, our models show perfect agreement with 
MNW13, indicating that the problematic trends at 2 = 
originate only at low-2. We conclude that compared to the 
old code version and to compilations of previous simulations 
(e.g lGuo et a"Lll2010l : ISawala et al.ll201ll . CS12, MNW13) our 
results look promising. 

We analyze the SFH problems further by separately 
plotting the SFHs of the most massive (M200 > 1.55 x 
IO^^Mq), least massive (M200 < 5. x IO^^Mq) and remain- 
ing haloes in Figure O We determine SFHs by following the 
in-situ growth of the main progenitor of the 2 = galaxy, so 
that accreted stars do not appear in the SFH. To be consis- 
tent with Figure [21 for the two 2 = major merger haloes, 
we add up SFHs of the two merging galaxies. In each of the 
panels we overplot the predicted average SFHs for typical 
haloes of the lowest, mean and highest 2 = halo mass 
in the panel as given in MNW13. These SFHs also do not 
include accreted stars. As discussed above, individual halo 
assembly histories can differ significantly from the assumed 
average assembly histories. The comparison can however re- 
veal trends in the SFHs of our sample of simulated galaxies. 
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Figure 4. Star formation rates SFR at 2 = 0, 2 = 1 and 2 = 
2 (averaged over 1 Gyr each) vs. stellar galactic mass Mstella r 
compared to the observed relations found by Elbaz ot alj 1120071 ) 
(z = 0, 1, black), [Paddi et al.l ||2007| ) (2 = 2, black) with Icr errors 
represented by dashe d lines and to median relations found by 
iKaiisawa et al.l 1I2OIOI ) (z=l,2, magenta and cyan). 



We start with the middle panel including AqC, which 
we have analyzed above. There is a trend for these haloes to 
overproduce stars at 2 > 4, with AqC showing the the high- 
est SFRs, which are however mostly caused by its early as- 
sembly of mass. Otherwise this trend agrees with the slightly 
high stellar masses at z ^ 2. As a population, these haloes 
show reasonable later agreement, leading to the agreement 
found in Figure [2] at 2^1. Probably due to the high-2 over- 
production, the predicted peak at 2 = 1 — 2 is not well repro- 
duced by the models. Moreover, there is a small tendency 
towards underproduction at 2 < 1. Generally, the conclu- 
sions for this sub-sample are the same as for AqC alone. 

For the low mass haloes in the left panel, there is hardly 
any overproduction at high-2, which might be connected to 
the lack of resolution and the consequently unresolved star 
formation in small progenitor haloes. Tests with higher reso- 
lution indicate that high-2 SFR increase with resolution for 
these haloes. The sub-sample agrees with the predictions at 
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Figure 5. Mean mass- weighted z = age of stars within 3 stellar 
half-mass radii i?5o vs. stellar gala ctic mass Mste ll ar com pared 
to the observed relation found by iGallazzi et al.l BOOST) for a 
magnitude-limited sample of galaxies (magenta lines). We also 
overplot in black th e relation for the lum inosity-weighted age of 
elliptical galaxies bv . Thomas et al.l 1120101) . Dashed lines represent 
1(T errors. 



2 ~ 0.5 — 5, but late SFRs are low. This leads to the agree- 
ment with abundance matching at 2 = 1, but to the low 
2 = masses seen in Figure [21 

The high-mass haloes in the right panel show a signif- 
icant overproduction of stars at high redshift only for halo 
AqA, which has a similar mass assembly history to AqC. 
There is good agreement with predictions at 2 ~ 0.8 — 4. 
However, SFRs are significantly too high at later times and 
lead to high 2 = stellar masses. We find that at these 
times there is still efficient infiow of gas onto the galaxies, 
but stellar feedback as applied in our models is too weak to 
efficiently remove gas from the disks as most of the ejected 
gas returns in galactic fountains. AGN, which we do not 
model, are generally believed to stop hot halo gas from 
cooling onto galaxies at these redshifts in massive haloes, 
though typically fo r higher masses than in our sample (e.g. 
ICroton et allbOOel ). 

To have a more direct comparison of our model SFRs 
to observation, we determine SFRs at 2 = 0, 1 and 2 for 
Figure [4] by averaging over 1 Gyr of star formation each. 
W e compare th e corre spon ding values to the re lations found 
bv lElbaz et"aLl (|2007l ) and iDaddi et al.1 (|2007l ) for observed 
SFRs vs. galaxy stellar masses. For the major merger mod- 
els 2283 and 0977 we here separately include the merging 
galaxies. Apart from our least massive halo 6782, in which 
there is hardly any star formation after 2 ^ 0.4, our model 
sample agrees well with observation and all model galaxies 
lie within 2cr of the observed relation at 2 = 0. The high- 
mass haloes do show a tendency for high SFRs, as expected, 
and the rest of the models show a slightly lower average 
than observed, as was also conclude d from Figu r e [3l I t has 
to be noted, that t he observations of lElbaz et al.l ((20071) and 
iDaddi et al.l (|2007|) only considered blue, star- forming galax- 
ies, whereas the results of MNW13 are an average over blue 
and red galaxies. 
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Figure 6. Edge- and face-on images of the z = stellar l^-band surface brightness Ev (left) and of the gas surface densities S for our 
model galaxies in haloes (top to bottom); 6782, 1646, AqB and one of the merging galaxies in 0977. 



At all three redshifts in consideration, the slope of the 
SFR vs. Msteiiar relation constituted by our models is in 
good agreement with observation. At z — 1, where the 
different observations agree well, SFRs are typically -^ la 
low, which we can connect to the fact that our models 
do not reproduce the peak in the SFHs at z — 1 — 2 
properly, but on average predict a flatter SFH. At z = 2, 
our model SFRs a r e '--^ 2a low compared to the results of 
iDaddi et al] (|2007D . lKaiisawa et all l |2010l ) find lower SFRs 
at Msteiiar < 10^^ Mq, which our lower mass galaxies agree 
with. The higher mass galaxies have significantly lower SFRs 
compared to these observations as well. This seems to dis- 
agree with the results from Figure O where no deficiency 
at that time is apparent. If we take into account that the 
model stellar masses at z = 2 are slightly high, it seems rea- 



sonable to assume the measured SFRs together with stel- 
lar masses according to the MNW13 predictions from abun- 
dance matching. This can explain part of the discrepancy, 
but it this is not enough to explain the full difference be- 
tween observations and models. Interestingly, iKannan et al.l 
(J2013h found a similar discrepancy with observations for 
their galaxies in a simulation of a cosmological volume. 

Another direct comparison to observation can be made 
for ages. We determine the mean age of populations within 
three stellar half-mass radii a nd compare them to t wo obser- 
vational relations in FigurelSl lGallazzi et al.l l|2005r ) provided 
an age-stellar mass relatio n for a magnitude-lim ited sample 
of SDSS galaxies, whereas Thoma s et al.l (J2010l ) examined a 
sample of morphologically selected elliptical galaxies. This 
figure very clearly reveals the trends discussed above and a 
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Figure 7. As in Figure[6]for galaxies in haloes (top to bottom): AqD, 0959, 0664 and 0616. 



division of our sample in three parts. The trend of ages of all 
model galaxies with stellar mass is almost flat and thus dis- 
agrees with the full galaxy population. The Aquarius haloes 
and the other haloes with masses ~ IO^^Mq are in rea- 
sonable agreement with both observational datasets, which 
strengthens the conclusions about their well-modeled SFHs. 
The massive galax;ies with high late SFRs only marginally 
agree with the youngest observed M* ~ 10^^ Mq galax- 
ies. The lower mass galaxy models, which include the three 
haloes with lowest masses and the galaxies undergoing a 
major merger at 2: = 0, overlap with the observed elliptical 
galaxies but are more than la high compared to the Gallazzi 
et al. sample, which at these masses is dominated by disks. 

In conclusion, our feedback model, calibrated on halo 
AqC, works well for all haloes with similar and slightly lower 
masses. These models show overly high SFRs only at high 



redshift. In the most massive haloes we studied, we find that 
SFRs are significantly too high at low redshift, possibly on 
account of the lack of a model for AGN feedback. For our 
lower mass galaxy models we find too low SFRs at 2: < 1, 
which makes the ages of their stellar populations too high 
for disk galaxies at these masses. 



5 MORPHOLOGY AND KINEMATICS 

After studying how the stellar mass assembles in our model 
galaxies, we now have a closer look at the morphology and 
the kinematics of their stellar components. 
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5.1 Structural Properties 

We start by a visu al inspection of gala x y mo rphologies. We 
use the models of iBruzual &: Charloj ( 20031 ) to assign V- 



band luminosities to the stellar particles according to their 
mass, age and metallicity. We then create edge-on and face- 
on surface brightness Ey images for the stellar component 
without taking obscuration into account. We also create gas 
mass images and present the results for eight model galaxies 
at 2; = in Figures [^ and [7] to display the variety of galaxy 
types produced in our simulations. 

In Figure [5] we present lower mass galaxies. All of these 
galaxies do s how extended gas disks, as do observed disk 
galaxies (e.g. IWalter et al.l 12008). Moreover, the gas disks 
and, with the exception of 6782, also the stellar components 
do show warps of varyin g extent. War ps are also frequently 
observed in real galaxies ( Sancisilll976l ') . Considering our full 
sample of 18 galaxies in the 16 simulated haloes (2 each for 
2283 and 0977), 10 galaxies show warps and 11 have gas 
disks that extend beyond the stellar populations. 

The stars of the galaxy in our least massive halo 6782 
(first row of Figure[6]) show an ellipsoidal distribution, which 
is only mildly flattened. The gas however lives in an ex- 
tended, warped disk, which is dominated by two spiral arms. 
The gas disk has formed after a major merger at z = 0.4. Our 
model 1646 (second row) also shows elliptical morphology 
in stars, the visible rotational flattening is however higher. 
Moreover, there are signs of a disk and a misaligned ring. 
The gas displays a strongly warped, relatively compact disk 
with disturbed spiral structure. We will show below that 
this galaxy features strongly misaligned infall at z ~ 0.1 
and a counter-rotating stellar population. Misaligned infall 
and reorientation of disks are frequent phenomena in our 
simulations. About half of the z = disks have experienced 
reorientation of 40 degrees or more since they started form- 
ing. 

Halo AqB (third row) hosts a very thin, extended stellar 
disk galaxy with flocculent spiral structure and a prominent 
bulge. The gas forms a very extended, slightly warped disk 
with spiral structure as in the stellar component. One of 
the disks in the major merger halo 0977 (fourth row) shows 
similarities to the one in AqB in terms of spiral structure 
and central bulge. However signs of the interaction with its 
neighbour are clearly visible. The system shows elongated 
tidal features and the stellar disk has been considerably 
thickened. The edge-on views reveal that the SF regions in 
the spiral arms do not lie in a well-defined plane as in AqB 
but in projection superpose to create a thick disk. 

In Figure [7] we display four galaxies with higher masses. 
0959 and 0616 show overly high SFRs at z = 0, whereas 
AqD and 0664 have reasonable z — masses. For these 
galaxies, only 0959 (second row) shows an extended, warped 
gas disk. It also displays a prominent stellar bar and a multi- 
arm spiral pattern. Only three of our 18 z = galaxies 
show clear signatures of bars, though more galaxies show bar 
signatures during their evolution . The z = rate of b ars is 
thus lower than in real galaxies (JEskridge et al.ll2000r ). The 
model galaxy in halo 0664 (third row) has a low gas fraction 
with gas being confined to a compact, un- warped disk The 
edge-on stellar morphology shows a disk and an underlying 
spheroidal distribution, the face-on image shows no distinct 




5x10" <M2oo/Ms< 1.55x10'' 




15 
R [kpo] 

Figure 8. Radially averaged z = face-on V-band surface bright- 
ness Sy as a function of radius R for all our model galaxies. We 
split the sample into three bins of halo mass: M200 < 5 X 10^^ Mq 
(top), 5 X 1O"M0 < M200 < 1-55 X IQI^Mq (middle) and 
M200 > 1.55 X IO^^Mq (bottom). 



spiral structure. The morphology is thus reminiscent of SO 
galaxies. 

In the next subsection, we will show that AqD has the 
highest kinematic disk fraction of all our model galaxies. 
The first row of Figure [7] reveals a large stellar disk with 
a small bulge and faint spiral structure. The gas lives in a 
thin, un-warped disk, which is less extended than its stel- 
lar counterpart. The morphology of the galaxy in halo 0616 
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(fourth row) is clearly also disky and features spiral arms. 
However, the disk appears significantly thicker and more 
compact. The gas disk is also compact and thick explaining 
the high SFRs. Moreover, the gas distribution extends ver- 
tically beyond the disk. We find it to be inflowing as part of 
a galactic fountain. 

Disk galaxies have long been known to show roughly 
exponential surface density profil es with central bulge con- 
centrations (|de VaucouleurdllQSa V In the recent years it has 
become clear that most of these exponentials actually consist 
of two or more parts which themselves are approximately ex- 
ponential. Most of these breaks (or truncations) in profiles 
are down-bending, however a signific ant fraction of galax- 
ies h as up-bending profiles (see e.g. iMartin-Navarro et al.l 
|2012| ). 

As the majority of our galaxies are disky their pro- 
files should show similar shapes, which indeed they do. In 
Figure [H] we plot radially averaged face-on "(/-band surface 
brightness profiles for all our model galaxies and split them 
into three panels according to halo mass (note the different 
scales). We find several interesting details: The lower mass 
galaxies are less extended tha n higher mass galaxies, as ob- 
served l|Courteau et al.l [20071 '). Also most galaxies do show 
(almost) exponential parts in their profiles. As in real disk 
galaxies, the central excess is more pronounced for higher 
mass galaxies and there are low mass galaxies showing no 
or hardly any central upturn. Four out of 18 galaxies have 
profiles that can be fit with a single Sersic profile with n ^ 1, 
one of which could be classified as 'pure-exponential'. 

Out of the galaxies with a central component, one has a 
single exponential outer component. In the remaining model 
galaxies, the disk components feature breaks. Nine of them 
are bulge-|-disk-f down-bending profiles, but four of our mas- 
sive galaxies (lower panel) show up-turning features in their 
outskirts as well. So from a qualitative point of view, the 
distribution of surface brightness profiles is in reasonable 
agreement with observations. 

Cosmological simulations of galaxy formation have long 
suffered fr om an over-concentrat ion of baryons in the centers 
of haloes (JNavarro fc Benall99lh . This is reflected by circu- 
lar velocity Kirc = ^GM{< r)/r curves that are strongly 
centrally peaked (see e.g. some models in CS12). Observed 
gas rotation curves of galaxies are h owever fairly flat or r ising 
outwards in low mass disk galaxies (|de Blok et al.ll2008l ). To 
reproduce this observation, it has been shown that feedback 
models have to selectively remove low angular momentum 
gas (e.g. iBrook et al.ll201ll ). 

In the top panel of Figure |9] we plot circular velocity 
curves for all our model galaxies as a function of the radius 
normalized to the stellar half mass radius R/Rso- Only the 
two galaxies in 0959 and 0616, which show a strong overpro- 
duction of stars at low z and feature compact, thick gas-rich 
disks show central peaks in Vdrc, but to a less severe extent 
than e.g. several models in CS12. The other galaxies indeed 
show very fiat circular velocity curves with detailed shapes 
depending on the details of their formation histories. In gen- 
eral, our feedback mechanism is thus capable of producing 
realistic central mass distributions. 

In the other panels of Figure [9] we show the fractions of 
the contributions of dark matter, stars and gas to the Vdrc 
curves. Clearly the centrally peaked models are dominated 
by stars. It is interesting to note that there is a contirmous 
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Figure 9. Spherically averaged z = circular velocity Vdrc = 
\/GMtoti< r)/r for all our models (top panel). The middle 
panels shows the fraction of the contributions to Vcirc by dark 
matter (DM) and stars Vdrc, dm/ Vcirc, stellar, where, Kirci = 



y'GMi{< r)/r. The bottom panel depicts 



s/K 



cirCjgas/ ''circ, stellar ■ 



transition in our galaxies from low-mass galaxies for which 
at all radii the gravitational potential is dominated by dark 
matter to more massive galaxies, which are star-dominated 
in the centers and for which dark matter becomes dominant 
at a radius that increases with galaxy mass. As far as gas 
is concerned, the baryonic contribution is dominated in the 
centers by stars for all galaxies. However for low-mass, more 
gas-rich galaxies, the gas contribution surpasses the stellar 
contribution at the outskirts of the galaxies. 

In summary, our model galaxies display a range of re- 
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Figure 10. Distributions of circularity e = Vrot/Vrot,max(S) for our model galaxies at 2 = 0. The panels for the mergers 2283 and 0977 
contain distributions for both galaxies. For 2283 the two galaxies show almost indistinguishable distributions. 



alistic morphologies in gas and stars. Unlike in many previ- 
ous simulations they are not too centrally concentrated and 
show no significant disagreement with observed galaxies in 
terms of radial surface density profiles and circular velocity 



5.2 Circularity distributions 

We now want to connect morphology to stellar kinematics. 
Circu larity e = Vl-ot/Kot,max(S) distributions (Abadi et al. 
[2003) are a tool that have been widely used to quantify stel- 
lar disks in simulations. Here Vrot is the rotational velocity 
of a stellar particle in a cylindrical coordinate system with 
the disk axis as symmetry axis. Particles on circular orbits 
thus show 6=1 and a thin disk results in a peak close to 
this value. A bulge will results in a peak around e = 0. We 



have computed such distributions for all of our models and 
present them in Figure [TOl 

We notice, that apart from 6782, the lowest mass halo, 
all our models show distinct disk peaks. 6782 shows an asym- 
metric distribution indicating some rotational support, but 
no stellar disk as we had observed in the first row of Figure 
[6l The spheroid results from a major merger at 2 ~ 0.4 and 
a lack of SF afterward. 4323 shows a small disk and a domi- 
nant bulge, whereas 4349 shows a broadened disk peak indi- 
cating a rather thick disk. Considering that all these haloes 
have overly low SFRs in the low-z era of disk formation, the 
lack of dominant thin disks is not surprising. 

For the haloes with z = major mergers, 2283 and 
0977, we plot distributions for both galaxies. In each case 
they show very broad disk peaks, indicating that these are 
mergers of disk galaxies with similar masses. The galaxies 
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Figure 11. Rotation to dispersion ratios Vrot/cz at z = as a function of the fraction of stellar galactic mass that has formed until the 
formation time iform of a star particle relative to the stellar galactic mass at z = 0. Our model galaxies are split into three bins in halo 
mass: A/200 < 5 X 10"Mq (left), 5 X 1O"M0 < M200 < 1.55 X 10^'^Mq (middle) and M200 > 1-55 X IQI^A/q (right). 



are still separated but the gravitational interaction has al- 
ready considerably heated the stellar disks. The thickening 
of one of the galaxies in 0977 is clearly visible in the fourth 
row of Figure [B] 

The most peculiar distribution of circularities is exhib- 
ited by 1646. It consists of two disks with opposite directions 
of rotation that live on top of each other. The co-rotating 
disk is more massive, but is also considerably thicker than 
the young and relatively thin counter-rotating disk. Addi- 
tionally, a distinct bulge peak is visible. This composition 
explains the disturbed and more elliptical morphology visi- 
ble in the second row of Figure [B] The model galaxy in halo 
1192 shows a disk peak with a wide tail to lower values of e, 
which is connected to a bar present throughout the evolution 
of the disk. 

Among the most massive haloes, 0664 has the lowest 
low-z SFR and a correspondingly small disk peak, which 
has formed after a merger a.t z ^ 0.7. The circularity distri- 
bution is dominated by a broad bulge peak resulting from 
the merger. This explains the SO like morphology seen in 
Figure [T] 0959, 0858 and 0616 show dominant disk peaks 
with tails to lower values of e indicating rather thick disks. 
The e distribution of 0959 is shaped by the presence of a 
bar (see Figure [T] whereas the disks in 0858 and 0616 are 
heated by minor mergers occurring throughout their evolu- 
tions. The bottom row of Figure [6] correspondingly shows a 
compact and thick stellar disk in 0616. 

The most dominant and thinnest disk peaks indicat- 
ing high fractions of thin, star forming disks can be found 
in the pan e ls for the Aquarius haloes. Of the haloes from 
lOser et al.l l|20ld ) only 1196 shows a comparably thin disk 
peak, which however only forms after a merger at z ~ 0.7. 
This is not surprising as the Aquarius haloes were selected 
to be prime candidates to host such galaxies. However, CS09 



were not able to completely confirm this conception in their 
study of all haloes, as their highest kinematic disk fractions 
were ~ 25%. Higher disk fractions were found for different 
codes for AqC in CS12, but the se galaxies did n ot compare 
well to observed disk galaxies. lOkamotd (|2013l ) presented 
the best Aquarius cosmological disk models so far for haloes 
AqC and AqD, their circularity distributions are however 
less disk dominated than ours. Considering the idealized, 
semi- cosmological models in AqA and AqC from AW13, our 
cosmological simulations produce disk fractions similar to 
the reference model ARef, but not as good as the best mod- 
els presented there. This is not surprising, as these models 
ignore the evolution at z > 1.3 when destructive mergers 
occur. 

In summary, our models show a large variety of kine- 
matic properties ranging from rotating spheroids, bulge- 
dominated disks and counter-rotating components over disks 
thickened by misaligned infall, bars and mergers to compact, 
massive disks and dominant thin disks. 



5.3 Disk Fractions 

Disk fractions in simulations are sometimes quanti- 
fied by mock p hotometric surface brightness decompo - 
sitions (see e.g. IXgertz et al.l I2OIII : iGuedes et al] 120111 ). 
IScannapieco et al.l (2010|) showed that this way of estima- 
tion yields significantly higher disk fractions than kinematic 
decompositions. In this work, we focus on kinematic disk 
fractions way as they constitute a more fundamental char- 
acteristic of the stellar population. 

To have a direct comparison to the models presented 
in CS12 we have calculated a distribution of circularities 
ev according to their definition based on the local circular 
velocity of a stellar particle. We then determined the fraction 
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Figure 12. The normalized number densities of the stars in six of our galaxy models at 2 = over the plane of circularity e vs. formation 
scale factor aformi "(^ i '*form/"max). The models are from top to bottom: 6782, 1646, AqB, one of the two merging galaxies in 0977, AqC 
and 0616. 



f, = ev > 0.8. For AqC we find /, = 0.55. All models in 
CS12 siiowed /e < 0.55, witli tlie liigliest values occurring 
in galaxies with centrally peaked rotation curves. Moreover, 
our ev distribution for AqC shows a higher ratio of disk 
peak to bulge peak than those galaxies. We thus conclude 
that our model produces a better disk than any of those 
models. 

In AW13 we argued that rotation-to-dispersion ratios 
are a good additional tool to characterize disks. We sort 
all stars in a galaxy as a function of age and bin them in 
equal-mass bins. We define a disk plane perpendicular to 
the total stellar angular momentum and calculate Kot as 
the mean tangential velocity and az as the rms velocity in 
the direction of the rotation axis for all stars in a certain 
age bin. In Figure [TT] we plot the ratio Kot/o"z as a function 
of the fraction of the z — stellar mass that has formed up 
to the time in consideration. We again split the models into 
three bins according to halo mass and include two galaxies 
per major merger halo. 

The idealized disk models from AW13 showed a mono- 
tonic decrease of Kot/o"z with increasing age. However these 
models did not include mergers, which for all haloes occur 
at high z. We thus expect an additional component with 
Viot/o'z ~ for the oldest ages. Moreover, the idealized mod- 
els showed that, at the resolution in consideration, values 



of Vrot /cTz » 10 as they are observ ed for the Milky Way 
thin disk (e.g. lAumer fc Binnevll2009h are prevented by the 
coarse resolution and the specifics of the modeling of ISM 
physics (see also lHouse et al.ll201lh . In general populations 
of stars with Kot /o"z > 2 can be interpreted as disky and 
the plot can thus be used to estimate disk fractions. 

Significant counter-rotating components are only found 
for the 15 % youngest stars in 1646 as discussed above and 
for the mildly (counter-) rot at ionally supported oldest bulge 
component of 4323. All models show dispersion-dominated 
bulge populations for their oldest components, as expected. 
For the low mass haloes only 4349 shows a significant 'thin' 
component with Kot/o"z > 5. The merger-thickened disks in 
haloes 2283 and 0977 show Kot/cz < 2 for all populations 
and the disk populations of 1192 and 1646, which suffer from 
a bar or misaligned infall, all show Kot/cTz < 5. So, as sug- 
gested in AW13, misaligned infall or reorientation of the 
disk plane heats existing stellar populations and prevents 
the formation of thin disks. 

All other haloes show distinct disk populations as dis- 
cussed for Figure [TOl The highest ratios Kot /o"z ~ 10 are 
displayed by the disks in AqB and AqD, but the massively 
star forming disks in the most massive haloes also reach 
Vrot/o'z ~ 8. If disks are present their Kot/o"z ratios de- 
crease almost monotonically with age, a sign for continuous 
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Figure 13. An overview of the disk fractions at z = in our model galaxies. In all panels we compare to our definition of thin disk 
fraction /thin given by the mass fraction in the 'thin disk area' of the e-a plane (see Fi£. ll2p . The left panel shows /{e > 0.7), the fraction 
of sta rs with circularity e > 0.7, the middle panel the fraction of the kinetic energy in rotation k = Erot / Ey^in a-s defined by I Sales et al.l 
1120121 ) and the right panel the scale factor when disk formation starts, Qdiskformation- 



disk heating (e.g. I Jenkins fc Binnevlll990l ). as discussed in 
AW13. The disks with the highest fractions of populations 
with Viot/o"z > 2 are AqA and AqD with fractions similar 
to 60 %. 

To gain a better understanding for the formation times 
of stellar disk components in our simulations we present in 
Figure [12] the distribution of stellar mass in the plane of 
circularity e vs. formation scale factor afoim- We present six 
different models to give an overview of how different evolu- 
tion histories show up in this plot. The SFH of 6782 is bursty 
and stops at z ^ 0.4, when a major merger occurs. Stars of 
all ages show a similar distribution of circularities, which 
characterizes an elliptical galaxy with small net rotation. 

A merger is still ongoing in halo 0977. Here we see two- 
stage formation. The stars that formed at z > 1.3 show 
a bulge distribution similar to 6782, whereas all younger 
stars live in a disk, which has been disturbed by the ongoing 
merger. The galaxy in 1646 shows three-stage formation. 
Again stars that formed at 2; > 1.3 live in a bulge, whereas 
stars that formed at 1.3 > z > 0.15 make up a thick disk. 
In this case the disk was heated by a drastic change in the 
angular momentum orientation of the infalling gas starting 
at 2; ~ 0.25. This reorientation results in the formation of a 
thinner counter-rotating disk. 

In 0616 the picture is more complicated. Stars at all 
formation times a < 0.8 show a broad bulge distribution. At 
a = 0.8 a last destructive merger occurs. However there is an 
older disk population, formed at 0.6 < a < 0.8, that has sur- 
vived thickened by the merger. A thinner disks forms after 
a = 0.8, but there are also young stars with thick disk cir- 
cularities, which are connected to ongoing interactions with 
satellite galaxies after a = 0.8. 

This is not the case for the disks in AqB and AqC. 
For both a dominant disk population is visible at high e 
after some well-defined time when formation starts, which 
for AqB is a destructive merger at 2 ~ 1. In AqC disk set- 
tling starts around 2 = 2 when mergers stop disturbing the 
galaxy. The disk populations widen in e for older ages, again 
indicating continuous disk heating. 

These e vs. a plots provide us with a good procedure 



to define disk fractions. We can define a 'disk area' in this 
diagram, which for AqB is clearly given by a > 0.55 and 
e > 0.8. The lower border for e is not well-defined, however 
e ~ 0.75 works well for all examples. We hereby also define a 
time when disk formation starts. With this definition, 6782 
and 0977-1 have no disks, the latter as the stellar system 
has been significantly heated by the ongoing merger. We 
therefore refer to this definition of disk fraction as 'thin disk 
fraction' /thin, but caution that observed thick disks would 
not necessarily be excluded by this method. 

For 1646, the counter-rotating disk is counted as the 
thin disk and for 0616 only the disk that formed at a > 0.8 
is included. We use this procedure to determine /thin for 
all model galaxies and use the values in Figure 1131 In the 
right panel we plot disk formation time aform against /thin- 
The values for /thin range between 0.15 (the counter-rotating 
disk in 1646) and 0.64 (AqA and AqD). We clearly see a 
correlation in the sense that younger disks have smaller mass 
fractions, which is not surprising. It is however interesting to 
have qualitative predictions for the connection of the time 
of the last destructive event and the final disk fraction. It is 
also interesting that the best disks in the Aquarius haloes 
AqA, AqC and AqD separate from the rest of the population 
indicating that their selection criteria indeed yielded haloes 
which host dominant disk galaxies. 

We also compare to previously used definitions to quan- 
tify disks. CS12 used all stars with ey > 0.8. We find that 
for our definition of the circularity e > 0.7 gives similar 
results. From the left panel in Figure [13] we see a clear cor- 
relation between the two disk fractions. However a simple e 
cut assigns a disk fraction to the t hick merging syste ms as 
well as to the elliptical galaxy 6782. ISales et al.l (|2012l ) used 
the fraction of kinetic energy in rotation k = -Erot/-Bkin as 
an indicator for diskiness. n increases with increasing /thin 
except for ongoing mergers. 

In conclusion, we have shown how rotation-to- 
dispersion ratios and the distribution of stars in the plane 
of circularity vs. formation time can be used to connect the 
kinematics and the formation history of model galaxies. Our 
galaxies are not as dynamically cold as real galaxies due to 
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numerical issues. However, they realistically show that the 
oldest components of galaxies are dispersion dominated and 
that rotational support increases with decreasing age of stars 
as observed in disk galaxies. We have shown that the model 
disk fractions in our Aquarius galaxies are higher than in 
any previously published simulation for these haloes. The 
disk fraction in our simulations depends mostly on the time 
since the last destructive event, be it merger or misaligned 
infall. 



6 SCALING RELATIONS 

We continue our study of the galaxy population constituted 
by our models and have a look at how they compare to 
observed relations involving stellar and gas masses, rotation 
velocities, sizes and metallicities. 



6.1 Gas fractions 



galaxies show a clear 
at lower gala xy masses 



Observed gas fractions in z - 

trend of higher gas frac tions 

iJHavnes fc Giovanellill2012l . see lPeeples fc Shank ar 201 1, for 
a rece nt compilation of data). At 2 = 1 — 2, IXacconi et al.l 
l|201[i ) using measurements of CO showed that gas fractions 
at a given stellar galactic mass are on average significantly 
higher. They however analyzed galaxies more massive than 
the ones in our simulations. Erb ct al. (2006) provided gas 
fractions for lower mass z ^ 2 galaxies by estimating the gas 
mass applying the local Kennicutt-Schmitt relation to SFR 
measurements. At higher masses, they foun d significantly 
lower gas fractions than those measured by iTacconi et al.l 
l|2010l ). 

These observations provide an important test for our 
models of star formation and feedback in the context of cos- 
mological gas accretion. Star formation uses up gas, whereas 
feedback can prevent gas from cooling and thus forming stars 
or eject gas from galaxies. Ejected gas can leave the halo if 
energetic enough or come back in galactic fountains. Clearly 
inappropriate modeling of these processes can thus be iden- 
tified from such a relation. 

In Figure [2] we attempt such a comparison. We deter- 
mine the mass of gas with T < 10* K within 20 kpc and 
use it to determine the fraction of cold gas to total baryonic 
mass in the galaxy. At z = 0, our simulations reproduce 
the observed decrease in gas fraction with increasing stellar 
mass which becomes shallow or potentially flat at masses 
M* > 10^°'^ Mq. For the shallow part of the relation, the 
gas fractions in our model galaxies scatter around ~ 15 %, 
as in real galaxies. For the lower mass galaxies, gas fractions 
can be as high as 65 % and the low-mass population in total 
shows a higher mean than observed. The low-mass galax- 
ies that lie significantly above the relation all live in haloes 
which show low stellar masses in the Mstdiar — M200 relation 
discussed in Figure [2] Thus for these galaxies the low SFRs 
at low-z lead to an underproduction of stars despite a signif- 
icant reservoir of gas in the model galaxies. This inefflcient 
SF then leads to overly high gas fractions at z = 0. 

As in real galaxies, the gas fractions in our model 
galaxies at fixed Mstoiiar steadily increase when going to 
redshifts 2 = 1 and 2. At 2 = 1, the gas fractions for 
Afstoiiar < W^^Mq are intermediate to the z — and z — 
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Figure 14. The fraction of the mass in gas with T < 10 K to 
total baryonic mass within 20 kpc vs. the galactic stellar mass 
A^stcllar a-t redshifts, z = 0, 1 and 2. Overplotted are the observe d 
mean z = relation presented by iPeeples &: Shankan 112011^ . 
an approx i mate relation for 2 = 2 obtained f r om th e results of 
lErb et al.l 1I2OO6I ') and data from lTacconi et al.1 1I2OIOI) . 



observational constraints. For Mgtdiar > 10^ ° M^, the gas 
fractio ns agree well with the observations of iTacconi et al.l 
l|2010r ). At 2 = 2, our models lie on top of the observational 
relation provided by lErb et al.l (|2006l ). with only a slight 
discrepancy towards low gas fractions for the lowest stellar 
masses, consistent with the fact that these galaxies show 
too high stellar masse s at these times. Con sidering the gas 
fractions measured by ITacconi et al.l l|2010l ) for higher mass 
galaxies, our ga s fraction in high e r mas s galaxies at 2 = 2 
appears too low. [Narayanan et al.l l|2012f) recently reported a 
similar discrepancy for the simulations of lDave et al.l l|2010l ) . 
They argue that this discrepancy arises from an overestima- 
tion of gas fractions based on incorrect CO-to-H2 conversion 
factors. 

In summary, the evolution of gas fractions in our model 
galaxies reproduce observed trends well. Deviations from ob- 
servational data can be linked to deviations from the pre- 
dicted SFH and uncertainties in high-2 observations. 
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Figure 15. Stellar half-mass radius R50 vs. galactic stellar mass 
Mgteiiar s-t redsliifts z = 0, 1 and 2. Overplotted in all panels is 
the ob served relation for late-type SDSS galaxies by IShen et alj 
1120031) for 2 = with Icr scatter denoted by dashed lines. 
iNagy et al.l 1 12011^ find that galaxies at 2: = 1.5 — 2.0 are typi- 
cally 27 % and galaxies at z = 2.0 — 2.5 typically 42 % smaller 
than galaxies of the same mass at 2 = 0. We overplot the corre- 
sponding lines in the z = 2 panel. 



6.2 Sizes 

Another important test for galaxy formation models is the 
size distribution of galaxies. Galaxy sizes depend on the dis- 
tribution of baryonic angular momentum and reproducing 
this distribution is crucial for achieving correct morpholo- 
gies, and thus also SFRs, gas fractions etc. Cosmological 
simulations have long suffered from an over-concentration 
of lo w angular momentum material in the centers of haloes 
([Nayarro & Benz 1991). It has been shown that in order to 
reproduce observed galaxy sizes the invoked feedback mech- 
anism has t o preferentially remove gas w ith low angular 
momentum ()Dutton fc van den Boschll2009l ) . Thus feedback 
modeling is effectively tested by galaxy sizes. 

In Figure [15] we attempt a comparison of o ur galaxy 
model sizes with observations bv lShen et al.l l|2010l ) for SDSS 
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Figure 16. Baryonic mass (within r = 20 kpc) vs. circular ve- 
locity Vcirc = \/GM(< r)/r at r = 4R50 at z ==0. Overplot- 
ted are the observational results from McGau ghl |2012h (black), 
Hall ct al. (2012) (magenta) and Avila-Rccsc ct al.l | |20()1 ') (cyan), 
for the so-called baryonic TuUy-Fisher relation, with la deviations 
represented by the dashed lines. 



late-type galaxies and bv lNaev et al.l IJ201 j l for z = 1.5 — 3.0 
star-forming galaxies. The half-mass radius -R50 of observed 
galaxies at all reds hifts increases on average with galaxy 
stellar mass Afstciiar.|Nagx.et_alJ (J2011h find that galaxy sizes 
at a given Mstdiar increase with time, being ~ 41% smaller 
than local galaxies a,t z — 2.0 — 2.5 and ~ 27% smaller at 
z = 1.5 — 2.0. This suggests that arou nd z = 1 galaxies start 
to lie on the local relation (see also Harden et al.ll2005l ). 

For our model galaxies, we also detect such an increase 
of R50 with time. At z = 2 the sizes scatter around values 
similar to those found by Nagy at z = 1.5 — 2.0, whereas 
at z = 1 and at 2: = they scatter around the z = rela- 
tion. At 2 = 1 — 2 the increase of R50 with stellar mass also 
agrees well with observations. At 2 = our models lie on an 
overly flat relation. We already showed that the most mas- 
sive haloes in our sample at z = host compact, overly star- 
forming galaxies. Here these galaxies explicitly appear as too 
small. On the other hand, our sample also includes two ma- 
jor mergers. The two galaxies in halo 0977 show elongated 
features due to their gravitational interactions and thus are 
most extended relative to the observations. Ignoring these 
model galaxies yields a relation with reduced scatter and a 
clear trend of increasing size with increasing mass, which is, 
however, still shallower than observed. 

Considering the connection of sizes and feedback mod- 
elling we note that our prescriptions are capable of repro- 
ducing observed galaxy sizes well at z ^ 1, but tend to be 
too strong at late times in removing low angular momentum 
gas in haloes with M200 < 0.7 x lO^^Af0 and too weak in 
haloes with M200 > 1-6 x IO^^Mq. 



6.3 The baryonic Tully-Fisher relation 

The luminosities and rotation velocities of disk galaxies are 
strongly correlated, they lie on the Tully-Fisher relation 
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||Tu11v fc Fisherlll972l '). In recent years, it has been shown 
that using the total baryo nic mass instead of luminosity 
yields a tighter correlation l|McGaughll2012l 'l. As McGaugh 
gives a relation that is among the steepest and tightest found 
in the literature, we also c onsid er the observatio nal results 
of lAvila-Reese et alj (|2008) and ' Hall et "al] (|2012l '). who find 
shallower slopes and significantly more scatter. We attempt 
to compare to these observational relations in Figure [TBI We 
plot the mass of stars and gas within r^^i (as defined in Sec- 
tion U21) against the circular velocity Vcirc = \/GM(< r)/r 
at r = 4 Rso- We use this radius, as iMcGaughl (|2012l ') use 
measurements of velocities where rotation curves are flat, 
which is well reproduced by our choice as depicted in Figure 
[5] Note that lAvila-Roosc ot al. (2008) use the maximum ro- 
tation velocity and Hall ct al. (20i2) use the velocity at half 
light radius. Only for the most massive, compact galaxies in 
our sample which overproduce stars at low z, would a dif- 
ferent choice of radius alter the data-points significantly. A 
different choice for the radius within which we consider the 
baryonic mass does also not change our conclusions. 

The baryonic TuUy-Fisher relation is intrinsically con- 
nected to the M* — Mhaio relation discussed in Figure [2] 
and to the gas fractions discussed in Figure [T4l Moreover as 
was shown e.g. in CS12, the correct distribution of mat- 
ter and thus the sizes of galaxies as discussed in Figure 
1151 also play an important role. The TuUy-Fisher relation 
is thus a good test for the interplay of the model ingredi- 
ents included in our galaxy formation code. Considering that 
for all the mentioned Figures we found at least reasonable 
agreement with observation, with smaller deviati ons we ex- 
pect a lso to find this in Figure 1161 Compared to iMcGaughl 
lJ2012f ). this is not, however, the case. Only a minority of 
galaxies lies within la of these observations and the power- 
law slope con stituted by our model sam ple is closer to 3 as 
ad vocated by lAvila- Reese et al.l (120081) and not 4 as found 
bv IMcGaughl (120121') . Due to the greater scatter found by 
iHall et al.l (|2012r ). our galaxies are marginally consistent 
with th eir results and the r e is g ood agreement with the re- 
sults of lAvila-Reese et all (|2008i l. 

Considering that the power-law slope of the relation 
given by our model galaxies is smaller than most obser- 
vations, it is interesting to discuss the trends. The most 
massive galaxies have too high rotation velocities for their 
masses. This agrees with the fact that they are too compact 
for their baryonic mass. There is a weaker trend for low- 
mass galaxies to lie above the relation, possibly connected 
to too efficient feedback for these haloes preventing a further 
concentration of the halo. 

In summary, our model galaxy sample shows a power- 
law slope for the baryonic TuUy-Fisher relation, that agrees 
only with the lowest in the range of values found observation- 
ally. The deviations from the relation can be connected to 
the remaining fine-tuning problems for the feedback mech- 
anisms which have already been discussed in the previous 
sections. 



6.4 Metals 

Another way to test how well the complex interplay of gas 
accretion, star formation, galactic winds and gas recycling 
is modeled, is to consider the amount of metals stored in 
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Figure 17. Mean mass- weighted galactic stellar metallicity 
within r = 3 -R50 , \Z\ vs. stellar galaxy mass Meteiiar at z = 
compared to the observational results of JGallazzi et al.l l l2005r) 



the galactic stars and gas. Metals are produced in stars and 
returned to the ISM in stellar explosions or winds. Thus 
metal production and energetic feedback are intrinsically 
coupled. Consequently, galactic winds are an efficient way to 
get metals out of galaxies. However, the metals in outflows 
can diffuse into hot halo gas and re-accrete, alternatively, 
outflowing material can return in galactic fountains, provid- 
ing ways for ejected metals to re-enter the galaxy. Moreover, 
the metal content of galactic gas can be diluted by the ac- 
cretion of pristine gas. Basically all ingredients in our galaxy 
formation simulations have a direct influence on the galactic 
metallicities. 

Whereas gas metallicities probe cumulative effects over 
the whole evolution of galaxies, stellar metallicities contain 
probes from every stage in galaxy formation. In Figure lTTI we 
attempt a comparison of the mean stellar met allicities of our 
model galaxies to the observational data of ICallazzi et al.l 
(2005). Apart from the too-metal- rich massive galaxies, the 
sample agrees well with observations and reproduces the 
slope of total stellar metallicity \Z\ with stellar mass. 

The picture is different however for gas phase metallic- 
ities, which we explore in Figure [181 The observed relations 
are quite uncertain in their calibration, which is why differ- 
ent authors find deviations of up to 0.3 dex in the normal- 
ization of the relation of oxygen abundance 12 + log{0/ H) 
versus galaxy stellar mass. The slope also differs, but too 
a smaller amount. Typically, quoted la errors at 2: = 
are of the order 0.1-0.2 dex depending on mass and cal- 
ibration. _We_jnclud«_three_differen^_obs^^ rela- 
tions b vlTremonti et al.l (I2OO4I ') . I Andrews fc Martinil (|2013l ') 
and by iMaiolino et al.l (|2008l ) and compare our results to 
them. Our galaxies follow a significantly steeper relation 
with hig her mass galax i es (M ^. > 2 x 10^°Mb) agree- 
ing with IMaiolino et all (|2008l ') and iTremonti et all (|2004l 'l 
and the lower ma s s gal axies in better agreement with 
[Andrews fc Martinil (|201j ). but on average more metal-poor 
than all three observational relations. 

At higher redshifts, we do not find such a disagreement 
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with the observation al slope of the relation . When we com- 
pare to the results of lMaiolino et al.l (|2008l ) at z = 0.7 — 3.5, 
we find reasonable agreement between our models and ob- 
servations when taking into account the observational un- 
certainties. At z = 3.5, our model galaxies are slightly too 
metal rich, which considering the over-production of stars 
at 2 > 4 is not surprising. At 2; = 2 we find perfect agree- 




Figure 19. Mean luminosity-weighted galactic stellar oxygen to 
iron ratio [O/Fe] vs stellar galaxy mass at z = compared to 
the observational results of [Johansson et al.l 1120121 ) for elliptical 



ment with observations, whereas at 2; = 0.7 our models are 
slightly metal-deficient, probably due to the small deficit in 
SFRs around z = 1. The discrepancies at 2 = obviously 
origin at z < 0.7. 

At late times, four lower mass galaxies are metal defi- 
cient. They include two of the haloes with the lowest masses 
and the two galaxies involved in the major merger of halo 
0977. We have noted that these galaxies are too gas rich. 
From observed relations, we know, that higher gas fractions 
and lower metallicities are present in lower mass galaxies. As 
our lower mass galaxies are deficient in SF and too gas rich, 
late low metallicity infalling material can dilute the metal- 
licity and lead to even lower metallicities than the observed 
ones. The most massive galaxies tend to be too metal-rich, 
which, as discussed for stellar metallicities, can be connected 
to overly high \ow-z SFRs. 

It is well known that SNII and SNIa produce different 
fractions of elements and explode on different timescales. 
Because of that, fractions of O (SNII) to Fe (SNIa) test the 
validity of the modeling of metal production yields, SFHs 
and model SN time distributions. In Figure [191 we plot the 
mean luminosity-weighted element r atio \0/Fe\ in stars. We 
compare them to observations from [Johansson et al.l (|2012l ) 
for elliptical galaxies. Although the average [O/Fe] in our 
galaxies agrees with ellipticals, these values are too high, 
as disk galaxies typically show lower values of [O/Fe] due 
to the higher SFRs at late times. Compared to observed 
ellipticals, the increase in [O/Fe] with galaxy mass is much 
shallower in our model galaxy sample. Clearly, higher low- 
2 SFRs in lower mass galaxies would shift [O/Fe] to lower 
values and the opposite is true for the highest mass galaxies 
which overproduce stars at low redshifts. However, to shift 
the normalization, a different calibration of yields and SN 
distributions over time is necessary. As we noted in Section 
12.21 the metal production yields we apply were not calibrated 
in any way. 

The modeling of turbulent diffusion of metals also has 
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an influence on the metaUicity in galaxies. We find that us- 
ing higher normalizations for the diffusion coefficients yields 
higher m etallicities in galaxi es. Using the normalization pro- 
posed by iGreif et al.l (|2009r ) (see Section 12. 3p increases the 
metaUicity of model galaxies by ^ 0.2 — 0.4 dex and would 
lead to significantly worse results compared to observations. 
The reason for this is, that outflowing gas is metal-enriched 
and thus loses metals to the circumgalactic gas by diffusion. 
The halo gas can later accrete onto the galaxy fe hen et al.l 
[2OI0I). This effect is enhanced by stronger diffusion coeffi- 
cients. 

Metals in galaxies are not equally distributed. Most 
galaxies show gas oxygen metallicities that decrease out- 
wards, with gradient slopes in units of d ex per scale le ngth 
being independent of galaxy mass (|Zaritskv et al.lll994r ). As 
more massive galaxies are more extended, absolute gradients 
in dex/kpc are thus steeper for galaxies with lower masses 
and smaller sizes. Disk galaxies have stee per gradients tha n 
ellipticals and than interacting disks (iKewlev et al.ll2010l ). 
IZaritskv et al.l (|l994l 'l found gradients ranging from -0.23 to 
-1-0.02 dex/kpc for a sample of ~ 40 disk galaxies with a 
mean ~ —0.06 dex/kpc. Gradients are believed to arise from 
the inside-out formation of disk galaxies and the correspond- 
ing variation of SF and feedback efficiency with radius. 

In Figure !^ we plot the gradients which we find in the 
oxygen abundances in the gas disks of our model galaxies 
and plot them against galaxy mass. All our galaxies show 
clear negative gradients, but they are rather shallow rang- 
ing from -0.053 to -0.007 dex/kpc , being a factor of ~ 3 — 4 
shallower than observations. We do however reproduce the 
fact that absolute gradients are steeper for smaller galaxies. 
As far as our interacting galaxies in the merger haloes 2283 
and 0977 are concerned, the extended disks in halo 0977 in- 
deed have flatter gradients than the other galaxies of similar 
mass. The significantly more compact galaxies in halo 2283 
however show the steepest gradients in our sample. 



It is interesting that IPilkington et al.l (|2012l ) found 
steeper gradients in a set of cosmological simulations with 
AMR and SPH codes. Considering our modeling there are 
two ingredients which can make gradients become too flat. 
Turbulent diffusion of metals is modeled on the scale of the 
SPH smoothing kernel and thus is not independent of res- 
olution. In test simulations, we flnd that stronger diffusion 
produces flatter metaUicity gradients, as diffusion intrinsi- 
cally acts against metaUicity gradients. Moreover, our two 
phase enrichment model, which spreads half of the metals 
into the hot gas phase, intrinsically yields a relatively non- 
local enrichment of the ISM, which also weakens gradients. 

In summary, our model galaxies reproduce the observed 
relation of stellar metaUicity vs. stellar mass well. The evo- 
lution of the gas phase oxygen abundance is also reproduced 
reasonably well at z J5 1. At low-z the problems with SFHs 
discussed in Section 321 lead to deviations from observation. 
The oxygen-to-iron fractions [0/Fe] in our models are too 
high and the metaUicity gradients in our gas disks are too 
shallow. This indicates that the details of our modeling of 
metal production, enrichment and diffusion still needs fine- 
tuning. 





0.01 r 




0.00- 







C) 


: 


V 


- 


X 


-0.01 r 


a? 




T3 


: 


r 


-0.02 r 


0) 


: 


-n 


- 


tc 


-0.03 - 


D) 


: 


C/l 


: 


CD 


-0.04 r 


I 


: 





: 




-O.Ob r 




-0.06 ^ 




10 



L z=0 






i 


r 


% 


)K ^ 


^ ^^ ^ 1 


)K 






^ ; 


; )K 






•\ 


)K 




)K 


\ 


; \ 


, , 1 




\ 



10 
M 



10 

stellar [l^d 



10^^ 



Figure 20. The (lack of) gas phase metaUicity gradients in the 
galactic disks at z = 0. Mean gradients in [O/ H] are plotted 
against galactic stellar mass Af stellar- 



T CONCLUSIONS 

We have presented an extensive update t o the mul ti phase 
SPH galaxy formation code by iScarmap ieco et al.l (|2005l , 
[200a). This includes updated metal cooling rates based on 
individual element abundances, a continuous modeling of 
metal production in SNIa and AGB stars, an increase in the 
star formation threshold and a model for the turbulent dif- 
fusion of gas phase metals. Moreover, we now use feedback 
prescriptions which spread energy in kinetic and thermal 
form and apply a model for the effect of radiation pressure 
from massive, young stars on the surrounding ISM. 

We calibrate our feedback model to match the SFH of 
the previously well studied halo Aquarius C (CS09, CS12) 
predicted by results from abundance matching techniques 
o ut to z = 4: (MNW13) . By doing so, we confirm the findings 
of IStinson et al.l ([2013) that modeling feedback from young 
stars can play a very important role in shifting star forma- 
tion to later times, in increasing the efficiency of feedback 
at high z and thus in improving cosmological simulations 
of galaxy formation. We find that in order to reproduce the 
inferred SFH, the infiuence of radiation pressure from young 
stars on the ISM has to be stronger at higher redshifts. We 
achieve this by invoking a model in which the higher gas 
velocity dispersions at early formation stages lead to larger 
star formation regions and thus to greater column densities 
of dust, which enhance the effect of radiation pressure. Our 
model still predicts an overproduction of stars at z > 4, 
however, pointing to inaccurate modeling of SF at early for- 
mation stages. 

It was noted in CS12 that there is a priori no reason why 
this halo should not be an outlier in terms of its SFH. How- 
ever, considering the many previous galaxy formation mod- 
els for halo AqC from CS12, our model not only improves 
the shape of the SFH, but also yields the most dominant 
disk galaxy in this halo known to us. The galaxy properties 
show good agreement with a variety of observations. 

The efficiency of the coupling between the radiation 
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field and the dusty gas is uncertain and some have argued 
that the approach of iHopkins et alj (|2011f ). which is similar 
to ou r modeling, overestimates the e ffects of radiation pres- 
sure iJKrumholz fc ThompsonI 120131 ) . The influence of our 
radiation pressure model on the SFH is, however, a key in- 
gredient for improving the agreement of simulated galaxies 
with observations. We conclude that any kind of feedback 
process acting before the explosion of SNe with effects sim- 
ilar to those in our model is likely to achieve this effect (see 
also IStinson et all |2013| : iBrook et al] [2012.) . 

We have applied our code to a sample of 16 haloes 
with virial masses M200 ranging from 1.7 x IO^^Mq to 
2.4 X 10^^ A/0, including h aloes which were sele cted to have 
no recent major mergers (Springel e t al.l I2OO8I ) and haloes 
with low-z mergers (Oser ct al. 2010|). As a population 
the resulting galaxies reproduce the observed M* — Mhaio- 
relation well at all redshifts z — — 4. There are however 
trends of too high stellar masses at high z, of low- mass haloes 
under-producing stars at low z, and of high-mass haloes 
overproducing stars at late times. 

Our model galaxies also reproduce observations of gas 
fractions, stellar half-mass radii, star formation rates and 
gas phase metallicities at 2: = — 3 reasonably well. At 
z — 0, our sample also shows agreement with observed stel- 
lar metallicities and the galaxies provide a reasonable fit 
to the baryonic TuUy-Fisher relation. Deviations from the 
observed scaling relations can mostly be connected to the 
problems with SFHs mentioned above. Considering these re- 
sults, one could argue that our finding that reproducing the 
SFH helps in reproducing a range of observations to some 
degree validates the abundance matching methods used by 
MNW13. 

The problems at the high mass end could well be ex- 
plained by the lack of a model for AGN feedback, which 
could prevent the cooling of large mas ses of halo gas at low 
z, wh ich occurs in our models (see e.g. lPuchwein fc Springell 
I2OI3I ). The low-mass haloes indicate that our model is far 
from perfect, as feedback in these galaxies is too efficient at 
preventing star formation at late times. Again, the question 
arises whether our calibration method is appropriate. As the 
peak in conversion efhciency M*/Mhaio appears to happen 
at slightly lower masses than that of AqC, it is justified to 
note that l ower halo ma s ses co uld be more appropriate for 
calibration. iBrook et al.l (|2012l ) argue that calibrating their 
galaxy formation model at a halo mass M200 < W^^'^Mq 
enables them to produce galaxies that reproduce a range of 
observations up to the peak in conversion efficiency. 

Compared to simulations with the previous code version 
(CS09) the resulting morphologies are significantly more 
disk-dominated. 15 out of 16 models show stellar disk com- 
ponents. Morphologies include spheroidal components, ex- 
tended gas disks, disks that are disturbed and thickened 
by mergers, warped gas and stellar disks, disks suffering 
from reorientation and misaligned infall, barred disks, a 
galaxy with counter-rotating disks living on top of each 
other, thin extended disks with varying bulge fractions 
and massive disks, which are compact and thick. Radial 
surface brightness profiles of our galaxies show pure ex- 
ponentials and bulge -I- exponential disk profiles, as well 
as up - and down-ben ding break s, as found in observa- 
tions IjMartin-Navarro et al.l l2012h . Moreover, apart form 
two massive galaxies, all circular velocity curves are flat and 



none show strong central peaks of the kind encountered in 
many previous simulations (e.g. some models in CS12). 

We show how distributions of stars over the circular- 
ity e vs. formation time Oform plane can be used to visual- 
ize formation histories and extract disk populations. This 
tool nicely reveals mergers and misaligned infall of gas. In- 
terestingly, if misaligned infall and/or the reorientation of 
existing disks occur, the signatures of these processes which 
heat stellar populations, are very similar to those observed 
for idealized, semi-cosmological models studied by AW13. 
These processes are of interest in the context of the for- 
mation of counter-rotating compo nents and 2a signatu res 
found in several elliptical galaxies (|Krainovic et al.ll2008n . 



Most of the disks in our model sample feature a well- 
defined start for disk formation, defined by a merger or a 
misaligned infall event. The later the last destructive event 
occurs, the lower the disk mass. These 'disk ages' all cor- 
respond to z < 2. It is interesting that the age inferred for 
the thin d isk population in the so lar neighbourhood is even 
older fsee lAumer fc BinnevI 12003 and references therein). 

The disk fractions derived from the e — aform plots range 
from 15 to 65 %. Some galaxies are thus disk dominated. 
However, there are real galaxies with even higher disk frac- 
tions. Considering the shortcomings of our SFHs at z > 4, 
and assuming disk formation times are correctly modeled, 
we can estimate disk fractions for our haloes under the as- 
sumption that SFHs follow the MNW13 predictions. For our 
three most disk-dominated galaxies, these estimations yield 
an increase in disk fraction from 59-65 % to 70-77 %. We 
thus conclude that high disk fractions are not problematic 
for ACDM haloes with sufficiently quiet merger histories. 
The Aquarius haloes show the highest disk fractions and 
the quietest merger histories, as was expected from their se- 
lection. However, previous simulations (CS09, CS12) failed 
to validate these predictions. 

Some of the problems of our simulations are connected 
to metallicities. Although our model galaxies reproduce the 
stellar metallicity - stellar mass relation, the z — gas phase 
metallicity - stellar mass relation is too steep, our alpha- 
element abundances are too high and the metallicity gradi- 
ents in our disks are too shallow. No optimization of yields 
has been attempted for our models so far, leaving room for 
improvement. Shallow metallicity gradients are likely to be 
connected to the details of our diffusion model and our two- 
phase enrichment scheme, which both could require modifi- 
cation. Another problem is that tests indicate unsatisfactory 
convergence of our results when increasing the resolution. 
Such problems are common for galaxy formation codes (see 
e.g. the discussion in CS12). One must take into account that 
important ingredients are 'sub-grid' and thus that a depen- 
dence on the lowest resolved length-scale is expected when 
resolution is increased without adjusting model parameters. 

In general, our simulations provide us with an interest- 
ing sample of simulated galaxies with realistic properties at 
least for MW-like halo masses M200 ~ 10^^ M©. The simula- 
tions predict many more observables that can be compared 
to the real galaxy population than are touched in this paper. 
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